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Abstract 

A  component-based  approach  is  introduced  for  fast  and  flexible  solution  of  parameter-dependent 
eigenproblems.  We  consider  a  shifted  version  of  the  eigenproblem  where  the  left  hand  side  operator 
corresponds  to  an  equilibrium  between  the  stiffness  operator  and  a  weighted  mass  operator.  This  permits 
to  apply  the  Static  Condensation  Reduced  Basis  Element  method,  a  domain  synthesis  approach  with 
reduced  basis  approximation  at  the  intradomain  level.  We  provide  eigenvalue  a  posteriori  error  estimators 
and  we  present  various  numerical  results  of  modal  analysis  of  structures.  We  are  able  to  obtain  several 
orders  of  magnitude  speed-up  compared  to  a  classical  Finite  Element  Method. 
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1.  Introduction 

In  structural  analysis,  eigenvalue  computation  is  necessary  to  find  the  periods  at  which  a  structure 
will  naturally  resonate.  This  is  especially  important  for  instance  in  building  engineering,  to  make  sure 
that  a  building’s  natural  frequency  does  not  match  the  frequency  of  expected  earthquakes.  In  the  case 
of  resonance,  a  building  can  endure  large  deformations  and  important  structural  damage,  and  possibly 
collapse.  The  same  considerations  apply  to  automobile  and  truck  frames,  where  it  is  important  to  avoid 
resonance  with  the  engine  frequencies.  Eigenproblems  also  appear  when  considering  wind  loads,  rotating 
machinery,  aerospace  structures;  in  some  cases  it  is  also  desirable  to  design  a  structure  for  resonance, 
like  certain  microelectromechanical  systems. 

With  improvement  in  computer  architecture  and  algorithmic  methods,  it  is  now  possible  to  tackle 
large-scale  eigenvalue  problems  with  millions  of  degrees  of  freedom;  however  the  computations  are  still 
heavy  enough  to  preclude  usage  in  a  many-query  context,  such  as  interactive  design  of  a  parameter- 
dependent  system.  In  this  paper,  we  present  an  approach  for  fast  solution  of  eigenproblems  on  large 
systems  that  present  a  component-based  structure  -  such  as  building  structures. 

For  the  numerical  solutions  of  partial  differential  equations  (PDE)  in  component-based  systems,  sev¬ 
eral  computational  methods  have  been  introduced  to  take  advantage  of  the  component-based  structure. 
The  main  idea  of  these  methods  is  to  perform  domain  decomposition,  and  to  use  a  common  model  order 
reduction  for  each  family  of  similar  components.  The  first  and  classical  approach  is  the  component 
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mode  synthesis  (CMS)  [1,  2],  which  combines  static  condensation  at  the  interdomain  level  with  eigen- 
modal  truncation  at  the  intradomain  level.  For  parameter-dependent  equations,  an  alternative  model 
order  reduction  is  to  use  reduced  basis  (RB)  approximations.  This  idea  was  first  introduced  in  the  re¬ 
duced  basis  element  (RBE)  method  [3],  with  a  Mortar  approach  for  domain  decomposition.  In  this  case 
the  RB  approximations  provide  a  model  order  reduction  with  certified  accuracy  over  a  specified  range 
of  parameters,  and  the  same  RB  spaces  can  be  used  for  similar  components  with  different  parameter 
values. 

In  [4],  a  static  condensation  RBE  (SCRBE)  approach  is  developed  for  elliptic  problems.  It  brings 
together  ideas  of  CMS  and  RBE  by  considering  standard  static  condensation  at  the  interdomain  level  and 
then  RB  approximation  at  the  intradomain  level.  In  an  Offline  stage  performed  once,  the  RB  space  for  a 
particular  component  is  designed  to  reflect  all  possible  function  variations  on  the  component  interfaces 
(which  we  shall  denote  ports);  components  are  thus  completely  interchangeable  and  interoperable.  During 
the  Online  stage,  any  system  can  be  assembled  from  multiple  instantiations  of  components  from  a 
predefined  library;  we  can  then  compute  the  system  solution  for  different  values  of  the  component 
parameters  in  a  prescribed  parameter  domain.  The  Online  stage  of  the  SCRBE  is  much  more  flexible 
than  both  the  Online  stage  for  the  standard  RB  method,  in  which  the  system  is  already  assembled  and 
only  parametric  variations  are  permitted,  and  the  Online  stage  of  the  classical  (non-static-condensation) 
RBE  method,  in  which  the  RB  intradomain  spaces  already  reflect  anticipated  connectivity. 

In  this  paper,  we  present  an  extension  of  the  SCRBE  to  eigenproblems.  The  new  aspects  are  the 
following.  First,  the  SCRBE  normally  takes  advantage  of  linearity,  which  is  lost  when  considering 
eigenproblems.  Hence  we  begin  by  reformulating  the  eigenproblem  using  a  shift  a  of  the  spectrum  in 
order  to  recover  a  linear  problem.  Finding  the  eigenvalues  is  then  performed  at  a  higher  level:  using  a 
direct  search  method,  we  find  the  values  of  the  shift  a  that  correspond  to  singular  systems.  Second,  we 
provide  a  posteriori  error  estimators  of  the  eigenvalues,  not  only  with  respect  to  RB  approximations  but 
also  in  the  context  of  port  reduction. 

In  the  context  of  CMS  approaches  for  eigenproblems,  out  method  provides  some  important  features: 
treatment  of  parameter-dependent  systems  (as  explained  above),  optimal  convergence,  and  port  reduc¬ 
tion.  The  classical  CMS  only  achieves  a  polynomial  convergence  rate  [5,  6]  with  respect  to  the  number 
of  eigenmodes  used  at  the  intradomain  level.  This  can  be  improved  to  an  infinite  convergence  rate  by 
using  overlaping  components  [6],  but  at  the  expense  of  losing  simplicity  and  flexibility  of  component 
connections.  Our  method  somehow  provides  an  optimal  trade-off  since  it  retains  the  interface  treat¬ 
ment  of  classical  CMS  allowing  flexibility  of  component  connections  -  while  achieving  an  exponential 
convergence  rate  with  respect  to  the  size  of  RB  spaces  at  the  intradomain  level. 

We  also  provide  port  reduction  so  as  to  increase  even  more  the  speed  up.  Recent  CMS  contributions 
consider  several  port  economizations  (or  interface  reduction  strategies):  an  eigenmode  expansion  (with 
subsequent  truncation)  for  the  port  degrees  of  freedom  is  proposed  in  [5,  7];  an  adaptive  port  reduc¬ 
tion  procedure  based  on  a  posteriori  error  estimators  for  the  port  reduction  is  proposed  in  [8];  and  an 
alternative  port  reduction  approach,  with  a  different  bubble  function  approximation  space,  is  proposed 
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for  time-dependent  problems  in  [9].  We  can  not  directly  apply  CMS  port  reduction  concepts  in  the 
parameter-dependent  context,  as  the  chosen  port  modes  must  be  able  to  provide  a  good  representa¬ 
tion  of  the  solution  for  any  value  of  the  parameters.  In  this  paper,  we  adapt  to  parameter-dependent 
eigenproblems  a  port  approximation  and  a  posteriori  error  bound  framework  introduced  in  [10]  for 
parameter-dependent  linear  elliptic  problems. 

The  paper  proceeds  as  follows.  In  Section  2,  we  present  the  general  eigenproblem  and  its  shifted 
formulation;  we  then  describe  the  static  condensation  procedure.  In  Section  3,  we  add  reduced  basis 
approximations  and  develop  a  posteriori  error  estimators  for  the  eigenvalues  with  respect  to  the  corre¬ 
sponding  values  obtained  by  the  “truth”  static  condensation  of  Section  2.  In  section  4,  we  introduce  port 
reduction  and  provide  as  well  a  posteriori  error  estimators  for  the  eigenvalues.  Finally,  in  Section  5,  we 
present  numerical  results  on  bridge  structures  to  illustrate  the  computational  efficiency  of  the  approach. 

2.  Formulation 

2.1.  Problem  statement 

We  suppose  that  we  are  given  an  open  domain  12  C  d  =  1, 2  or  3,  with  boundary  <912.  We  then 
let  X  denote  the  Hilbert  space 

0}  , 

where  <912 £>  C  <9f2  is  the  portion  of  the  boundary  on  which  we  enforce  Dirichlet  boundary  conditions.  We 
suppose  that  X  is  endowed  with  an  inner  product  and  induced  norm  ||  •  ||  y-  Recall  that  for  any 

domain  O  in  Rd, 

Ft1(C2)  e{»£  L2{0) :  Xv  €  ( L2((D))d },  where  L2(0)  =  {z>  measurable  over  O:  f  v2  finite}. 

Jo 

Furthermore,  let  Y  =  L2(Xl). 

We  now  introduce  an  abstract  formulation  for  our  eigenvalue  problem.  For  any  p  £  V,  let  a(-,  •;  pi)  : 
X  x  X  — >  R,  and  m(-,-\p)  :  X  x  X  — >  R.  denote  continuous,  coercive,  symmetric  bilinear  form  with 
respect  to  X  and  Y,  respectively.  We  suppose  that  X^  C  X  is  a  finite  element  space  of  dimension  A f . 
Given  a  parameter  p  £  T>  C  Rp,  where  T>  is  our  parameter  domain  of  dimension  P,  we  find  the  set  of 
eigenvalues  and  eigenvectors  (A (p),u(p)),  where  A (p)  €  R>o  and  u(p)  €  satisfy 

a(u(p),v,p)  =  X(p)m(u(p),v;p),  Vu  G  X^ ',  (1) 

m(u(p),  u(p);  p)  =  1.  (2) 

We  assume  for  simplicity  that  the  eigenvalues  A n{p)  are  distinct,  of  mutiplicity  one,  and  sorted  such  that 

0  <  Ai  (p)  <  A 2  (a*)  ■  •  ■  <  A jv(/z). 

We  now  define  a  surrogate  eigenvalue  problem  that  will  be  convenient  for  subsequent  developments. 
For  a  given  “shift  factor”  a  €  R>o,  we  modify  (1),  (2)  such  that  for  any  p  €  22,  we  find  r(p,  u)eR  and 
\{p,  cr)  €  that  satisfy 

B(x{v,v),v-,p;<t)  =  T(p,a)a(x(n,<r),v-,p),  Vu  G  XM , 


a{x{P,°),x{P,a)\p)  =  1- 
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(3) 

(4) 


Here 


B(w,  v;  /x;  a )  =  a(w ,  v;  /x)  —  am(w,  v;  p) 


(5) 


is  our  “shifted”  bilinear  form.  Note  that  we  change  the  bilinear  form  on  the  right  hand  side  from  m(-,  •) 
to  a(-,  •),  which  corresponds  to  a  different  norm.  This  choice  is  motivated  by  error  estimation,  presented 
later  in  the  paper,  as  it  permits  to  derive  relative  error  estimates  for  the  eigenvalue  A„(/x). 

We  also  sort  the  set  of  eigenvalues  such  that  ti(/x,ct)  <  72  (/tx,  a) . . .  <  T^f(p,a)  -  note  that  due  to 
the  shift  the  first  eigenvalues  can  now  be  negative.  It  is  also  clear  that  %n( /x,  a)  =  ,  1  un{p)  for  any 

it  £  R,  so  we  shall  henceforth  write  y;n(/x).  Also 


so  that 


A n(p)  -  O 
A„(/x) 


(6) 


>0, 

if  0  <  a  <  A n(/x), 

(7) 

Tn(H,v) 

=  0, 

if  o'  =  An(/x), 

(8) 

<0, 

if  CT  >  Xn(p), 

(9) 

for  n  =  1, . . . ,  M . 


2.2.  Static  Condensation 

We  now  move  to  the  component  level.  We  suppose  that  the  system  domain  is  naturally  decomposable 
into  I  interconnected  parametrized  components.  Each  component  i  is  associated  with  a  subdomain  ft, , 
where 

i 

O  =  Hi,  f \  fl  flj/  =  0,  for  i^i!  . 

i- 1 

We  say  that  components  i  and  i'  are  connected  at  global  port  p  if  f \  n  Hj/  =  Tp  ^  0.  We  also  say  that 
7 1  =  Tp  and  jt,  =  Tp  are  local  ports  of  components  i  and  i'  respectively.  We  denote  by  nr  the  total 
number  of  global  ports  in  the  system,  and  we  denote  by  nj  the  total  number  of  local  ports  in  component 
i.  We  assume  that  the  FE  space  X1^  conforms  to  our  components  and  ports,  hence  we  can  define  the 
discrete  spaces  X^f  and  Z^j  that  are  simply  the  restrictions  of  X-Kf  to  component  i  and  global  port  p.  For 
given  i,  let  X{^0  denote  the  “component  bubble  space”  —  the  restriction  of  X ^  to  fli  with  homogeneous 
Dirichlet  boundary  conditions  on  each  7?,  1  <  j  <  nj , 

Xifo  =  {«|ni :  v  G  XM;  v\7j  =  0,  1  <  j  <  n]}. 

We  denote  by  A/J  the  dimension  of  the  port  space  associated  with  global  port  p,  and  we  say  that 

the  global  port  p  has  A f£  degrees  of  freedom  (dof).  For  each  component  i,  we  denote  by  k'  a  local 

port  dof  number,  and  I\i  the  total  numbers  of  dof  on  its  local  ports,  such  that  1  <  k'  <  K,  .  We  then 

introduce  the  map  Vi(k')  =  (p,k)  which  associate  a  local  port  dof  k!  in  component  i  to  its  global  port 

representation:  global  port  p  and  dof  k,  1  <  k  <  Afp  . 

To  formulate  our  static  condensation  procedure  we  must  first  introduce  the  basis  functions  for  the 

port  space  Z x  as  {£p,i,  •  •  •  ,  Cp,Nr}-  The  particular  choice  for  these  functions  is  not  important  for  now, 
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but  it  becomes  critical  when  dealing  with  port  reduction  we  refer  to  Section  4.  For  a  local  port  dof 
number  k!  such  that  Vi(k')  =  (p,k),  we  then  introduce  the  interface  function  £  XX ,  which  is  the 
harmonic  extension  of  the  associated  port  space  basis  function  into  the  interior  of  the  component  domain 
fl, .  and  satisfies 


•  V'c  =  0, 


Vu 


c  yW 

€  Ai;0> 


(10) 


V’fe' 


C p,k  on  Tp 

0  on  'yj  ^Tp,  1  <j<  nj. 


(11) 


If  components  i  and  j  are  connected,  then  for  each  matching  local  port  dofs  ki  and  kj  such  that  Vi(ki)  = 
Vj(kj)  =  (p,  k ),  we  define  the  global  interface  function  4>Pjfc  £  as 


ipk.  on 

P -k  =  \  ''Plj  on  -hl  (12) 

0  elsewhere  . 


We  will  now  develop  an  expression  for  Xn{p)  which  just  involves  dof  on  the  ports  by  virtue  of 
elimination  of  the  interior  dof  given  that  <j  =  A n(p)  -  starting  from  (13)  to  finally  arrive  at  (18).  Let 
us  suppose  that  we  set  a  =  Xn(p)  (for  some  n)  so  that  the  right-hand  side  of  (3)  vanishes.  Then,  for 
Xn{k)  G  X^  we  have 

B(Xn{p),  v;  p\  a)  =  0,  for  all  v  £  X1^ . 


We  then  express  XnilA  £  X^  in  terms  of  “interface”  and  “bubble”  contributions, 


/  nr  K 

Xnfa)  =E  bi(n,a)  +  EE  £lp,fc(/r,  crjUIp^,  (1^) 

i= 1  p— 1 fc=l 

where  the  UPtk(p,  cr)  are  interface  function  coefficients  corresponding  to  the  port  p,  and  bi(p,  cr)  £  Xf^. 
Here  is  independent  of  cr,  but  we  shall  see  shortly  that  we  will  need  and  Uk,P  to  be  cr-dependent  in 
general. 

We  then  restrict  to  a  single  component  i  to  obtain 


Bi(Xn(v),v,p.-,a)  =  0,  for  all  v  £  X^f0, 


(14) 


where  Bi(w,  v;  p\  a)  =  at{w ,  v\  p)  —  a mt(w,  v;  p),  and  where  a,;  and  to*  indicate  the  restrictions  of  a  and 
m  to  f 2j,  respectively.  Substitution  of  (13)  into  (14)  leads  to 

Ki 

Bi(bi{p,a),v;  p-,a)  +  ^2UVi{k)(p,a)Bi(ipitk,v,  p-,a)  =  0,  (15) 

k= 1 

for  all  v  £  Xff0. 

It  can  be  shown  from  linearity  of  the  above  equation  that  we  can  reconstruct  bi(p,  cr)  as 

Ki 

bi(p ,  cr)  =  ^2,Uv.{k){p,a)b^k{p,(j), 
k= 1 
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where  &i,fc(/ix,  cr)  €  Xff0  satisfies 


B,(bi,k(n,cr),v;n;<j)  =  -Bi(ipitk,v;  p;  a),  \/v  G  X{f0. 

Let  (A i,n(/i),Xi,n(At))  G  M  x  Xjf0  denote  an  eigenpair  associated  with  the  n  local  eigenproblem 
ai(Xi,n(tJ'),v;  H)  =  V v  G  xfi, 


then,  since 


.  Bi(v,v;fj,;a)  _  .  at(v,  v;  p)  -  crro;(u,  v;  p) 

»eX"  \Mh  v^C{f0  IMlx.i 


.  ai(v,  v;  p)  —  am.i{v,  v;  p)  .  rrii(v,v;p ) 

>  mi  - 7 - r -  mi  — - — rr^ - 


vex"  mi(v,v;n) 

/X  /..X  f  mi(v,v;n) 


;0  ll‘'NX,i 


=  (Xiti(p)-a)  inf 


vex"  ||f  iiA-,i 

the  bilinear  form  £%(■,  •;  /x;  er)  is  coercive  on  X^q  if  a  <  Xiti{p),  where  Aj)i(/i)  is  the  smallest  eigenvalue  of 
(17).  Hence  (16)  has  a  unique  solution  under  this  condition.  Note  that  we  expect  that  Xiti(p)  >  X\(p), 
and  even  Ajii(/z)  >  A n(p)  for  n  =  2  or  3  or  4  —  of  course  in  practice  the  balance  between  An  and  Aj>n/ 
will  depend  on  the  details  of  a  particular  problem. 

Now  for  1  <  k  <  Afp  and  each  p,  let 

®ptk(fJ,,cr)  =  Vp'k  +  ^  h,k'(p,v), 

i,k'  s.t.Pi(k')=(p,k) 

and  let  us  define  the  “skeleton”  space  Xs(p,cr)  as 

Xs(p,a)  =  span{$Pife(^,  a)  :  1  <  p  <  nr,  1  <  k  <  Up}. 

r 

This  space  is  of  dimension  nsc  =  1Up.  We  restrict  (13)  to  a  single  component  i  to  see  that  for 
a  =  Xn(p)  we  obtain 


x„(e)k  =  E  uVi(k)(n,  <*)  {h,k{P^a)  +  V’i.fe)  ■ 


This  then  implies 


x»(/‘)  =  EE  Up,k{n,cr)  §P,k{p,  cr)  S  Xs(p,a).  (18) 

p=i fc= l 

Then,  for  a  =  Xn(p)  and  /i  G  P,  we  are  able  to  solve  for  the  coefficients  UPtk{p,  cr)  from  the  static 
condensation  eigenvalue  problem  on  Xs{p,a):  find  Xn(p)  G  Xg(p,a),  such  that 

B{Xn(p),w,p;<r)  =  0,  Vu  G  X$(p,  a),  (19) 

a(Xn{p),Xn(p);p)  =  1.  (20) 

We  now  relax  the  condition  a  =  Xn(p)  to  obtain  the  following  problem:  For  a  G  [0,<Tmax]  and  p  G  V, 
find  (r„(/i,<r),xn(/b<r))  €  (M,  Xs{p,<j)),  such  that 

<3(xn(p,cr),v;p;a)  =  rn(p,  <r)a(xn(p\  cr),  v;  p),  \/v  e  Xs{p,a),  (21) 

a(xn(MW),Xn(MW);M)  =  1-  (22) 


It  is  important  to  note  that  this  new  eigenproblem  (21)  (22)  differs  from  (3)  (4)  in  two  ways:  first,  we 
consider  a  subspace  Xg(p,a)  of  X.j  and  as  a  consequence  t„(/x,ct)  >  Tn(p,a);  second,  the  subspace 
Xs{n,cr),  unlike  X-Kr .  depends  on  a,  and  furthermore  only  for  a  =  \n  does  the  subspace  Xg{p,a) 
reproduce  the  eigenfunction  Xn(p).  We  now  show 

Proposition  2.1.  Suppose  a  <  Ai,i(/x)  for  each  1  <  i  <  I  to  ensure  that  the  static  condensation  is 
well-posed. 

(i)  t„(/x,  a )  >  Tn(p,a),  n=  1, . . . ,  dim(Xs(fi,  a)), 

(ii)  t„(h,  a )  =  0  if  and  only  if  a  =  A„(/x), 

(Hi)  <j  =  A n(p)  if  and  only  if  there  exists  some  n'  <  n  such  that  rn>{p,  a)  =  0. 


Proof,  (i)  The  case  n  =  1  follows  from  the  Rayleigh  quotients 


Tl(^,cr) 


.  B(w,w;n;a) 

mf  — - - — , 

wex*  a(w,w\p) 


(23) 


and 


Ti(p,a) 


.  B(w,w;n;a) 

mf  — 7 - 7—, 

w£Xs(n.°)  a(w,w,p) 


(24) 


and  fact  that  Xs(n,  a)  C  X^ . 

For  n  >  1,  the  Courant-Fischer-Weyl  min-max  principle  [11]  states  that  for  an  arbitrary  n-dimensional 
subspace  of  X^r ,  Sn,  we  have 


ringer)  =  max 

W£Sn 


B(w ,  w;  /x;  a) 
a(w,  w,  p) 


>  Tn(p,a). 


(25) 


Let  Sn  =  span{xm(/x,  a),  m  =  1, . . . ,  n}  C  Xg{p,  a).  Then  r/n(p,  a)  =  Tn(p ,  a),  and  the  result  follows. 

(ii)  This  equivalence  is  due  to  (8). 

(iii)  (*t=)  Suppose  a  =  Xn(p)  for  some  n,  then  by  construction  %„(/x,  a)  £  Xg{p,a).  Since  the  same 
operator  B  appears  in  both  (19)  and  (21),  it  follows  that  Xn{p,cr)  is  also  eigenmode  for  (21),  (22)  with 
corresponding  eigenvalue  0.  That  is,  for  some  n' ,  Tn'(p,a)  =  0  is  an  eigenvalue  of  (21), (22).  Moreover 
Tn+i(p,  a)  >  Tn+i(p,  a)  >  0  so  n!  <  n.  We  note  that  if  n  >  1,  we  do  not  necessarily  capture  Xm{p)  m 
Xs(p,  a)  for  m  =  1,  2, . . . ,  n  —  1,  hence  it  is  possible  that  n!  <  n.  On  the  other  hand,  if  n  =  1,  then  we 
must  have  CT)  =  Xi(/b  a)- 

(=»)  Suppose  t n' (/x, (7 )  =  0  for  some  index  v! .  Then  X„/(yx,cr)  satisfies  (19),  (20),  or  equivalently, 
(3),  (4)  for  Tn (/x,  <7 )  =  0.  From  part  (ii)  of  this  Proposition,  this  implies  that  a  =  A„(/x).  Moreover 
t„+i(tx,  cr)  >  t„+i(tx,  cr)  >  0  so  n!  <  n. 

□ 


A  stronger  result  that  will  be  important  for  the  following  can  be  obtained  assuming 
Hypothesis  2.1.  For  1  <  n  <  nsc  the  functions  a  — >  Tn(p,a)  have  exactly  one  root. 
Corollary  2.1.  If  Hypothesis  2.1  holds,  <7  =  An(/x)  if  and  only  ifrn(p,a)  =  0. 
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Proof.  From  Proposition  2.1  (iii) ,  we  have  that  a  =  Ai(/z)  if  and  only  if  t\  (h,<t)  =  0.  Then  using 
Hypothesis  2.1  and  Proposition  2. 1  (iii) ,  it  is  easy  to  show  by  recurrence  that  a  =  A„(/z)  if  and  only  if 
r„(/z,  a)  =  0  for  1  <  n  <  nsc.  □ 


Lemma  2.1.  We  have  that 


for  each  n  =  1,2,... ,  A f . 


dTn{p,a)  _  1 

da  A 


Proof.  We  set  v  =  Xn{p)  in  (3),  to  obtain 


(26) 


13(Xn{p),Xn(fJ-)-,Pw)  =  r„{n,  a)a(Xn(tl),Xn(n)\tl)  =  Tn{p,  cr). 


(27) 


where  we  employed  (4)  in  the  last  equality  above.  We  note  that  Xn(ff)  is  independent  of  cr,  and  differ¬ 
entiate  with  respect  to  a  to  obtain 


d Tn{n,o)  d  .  . 

- -  =  ^%»(W,ln(W;/lA) 

d 

=  [a(Xn(p),Xn{tiy,  fj.)  -  am(xn(p),  Xn(p)-,  M)] 
=  -rn(xn(n),Xn(p);n) 

m{un(n),un{n)\n)  =  - 


Xn{t^) 


□ 


However,  the  result  does  not  apply  to  r„(/x,cr):  we  cannot  apply  the  argument  from  the  proposition 
to  (21),  (22)  since  in  general  y„(/i,  a)  depends  on  cr.  We  can  still  state  the  following 

Proposition  2.2.  Assuming  that  r„(/q  •)  is  differentiable  at  A n(n)  and  Hypothesis  2.1  holds ,  then 

drjjff  An)  _  1 

da  A n(p)' 


Proof.  We  know  that  r„(/x,  A „(/z))  =  Tn(p,  A „(/x))  =  0,  tPjfflPAl  =  and  rn(/x,  cr)  >  rn(/x,  cr).  So 


we  have 


Vh  <  0,  +  1 


V/i  >  0,  ^A»M  +  ft)>_ 


An(/x) 

1 


Since  Tn(p,  •)  is  differentiable  at  Xn(p),  we  have 

drn(n,  Xn(ii))  ..  rn(/z,  A„(/x)  +  ft)  .. 

- - -  =  lim  - - - =  lim 

da  h—>o~  h  h->-o+ 


A  n(/Li)’ 
rn(n,  An(/n)  +  h) 


Xn  (/^) 


□ 


To  assemble  an  algebraic  system  for  the  static  condensation  eigenproblem,  we  insert  (18)  into  (21), 
(22)  to  arrive  at 

nr  K 

EE  UPtk{v,  cr)S($p,fc(^,  cr),  u;  n;  cr) 

P= 1  k—1 

nr  Np 

=  r(^,cr)  EE  UPlk(n,  cr)o($p)fc(/u,  a),v ;  \i\  a),  Vu  G  Xs,  (28) 

P=  1  fc=l 
nr 

X]  E  UpAVi  Cr)a($p,fe(M,  O'),  $p,fc(/b  O');  M!  O')  =  1.  (29) 

p=l  fe  =  l 

We  now  define  our  local  stiffness  and  mass  matrices  A *(/u,  cr),  M*(/z,  cr)  G  for  component  i,  which 

have  entries 


A lk,k{p,,a)  =  ai(ipitk  +  bitk(v,o-),  ipi,k'  +  cr);  A*), 

Mk,ik(n,a)  =  rn (ipi, k  +  bi, k in.  O'),  4>i, w  +bi,k'{n,A:H), 

for  1  <  k,  k'  <  Ki.  We  may  then  assemble  the  global  system  with  matrices  B(/i,  cr),  A(/z,  cr)  G  Rn«xn«, 

r 

of  dimension  nsc  =  X)p=i-^p":  giyen  a  cr  G  K.  and  /i  G  T>,  we  consider  the  eigenproblem 

B(/x,a)Y(/x,(j)  =  T(n,a)A(n,a)Y(iJ,,a),  (30) 

Y(p,,a)TA(p,,a)Y(p,,a)  =  1,  (31) 

where 

M(n,a)  =  A(n,a)  —  (32) 

Assuming  (2.1)  holds,  we  can  recover  A „(/z)  as  the  value  of  a  for  which  the  nth  eigenvalue  of  (30)  is 
zero.  In  practice  we  use  Brent’s  method  [12],  which  is  a  combination  of  different  direct  search  methods 
for  root  finding  (bisection,  secant  and  inverse  quadratic  interpolation). 

3.  Reduced  Basis  Static  Condensation  System 

3.1.  Reduced  Basis  Bubble  Approximation 

In  the  static  condensation  reduced  basis  element  (SCRBE)  method  [4],  we  replace  the  FE  bubble 
functions  h;.fc(/qcr)  with  reduced  basis  approximations.  These  RB  approximations  are  significantly  less 
expensive  to  evaluate  (following  an  RB  offline  preprocessing  step)  than  the  original  FE  quantities,  and 
hence  the  computational  cost  associated  with  the  formation  of  the  (now  approximate)  static  condensation 
system  is  significantly  reduced.  We  thus  introduce  the  RB  bubble  function  approximations 

bitk(ix,a)  s=s  £>»,*,  (/x,  a)  (33) 

for  a  parameter  domain  (/z,cr)  G  V  x  [0,crmax],  where 

Umax  =  Co- min  min  A 
pex>i  <i<i 
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Here  eCT(<  1)  is  a  “safety  factor”  which  ensures  that  we  honor  the  condition  a  <  Aj)i(/x)  for  all  1  <  i  <  I. 
Next,  we  let 

^p,k(^a)  =  ^P,k  +  X!  Kki(n,(j), 

i,kiS.t.Pi(ki)=(p,k) 

and  define  our  RB  static  condensation  space  X$(/x,  a)  C  X-^  as 

Xs{h,(t)  =  span{$P;fc(^,  cr)  :  1  <p<  nr,l  <  k  <  A/J }. 

(Note  that  Xg(n,a)  qL  Xs(/j,,a).)  We  then  define  the  RB  eigenproblem:  given  (/n,  ct)  G  V  x  [0,<7max], 
find  the  eigenpairs  (r„(/x,  cr),  Vn(/z,  a))  that  satisfy 

I(/z,o-)¥(/i,o-)  =  r(/c,  a)A(/i,  a)V(/j,  a),  (34) 

V (/i,  a)TA(fx,  a)V (/j,,  a)  =  1,  (35) 

where  B(/z,  cr),  A(/z,  cr)  are  constructed  component-by-component  from 

A‘k,k([i,a)  =  +  bi,k(Li,o-),ipi,k'  +bi,k'(fJ-,cr);n), 

Mlk,tk(fi,<j)  =  mi(ipi}k +  b^k(n,  a),ipi}k'  +  bitk>(fi,,a);n), 

for  1  <  k,  k'  <  Ki:  and  where 

B4(^,  a)  =  A*(p,  cr)  —  crM*(/U,  a).  (36) 


3.2.  Reduced  Basis  Error  Estimator 

We  now  consider  error  estimation  for  our  RB  approximations.  First,  since  Xs(/z,  cr)  C  X^ ,  by  the 
same  argument  as  part  (i)  of  Proposition  2.1,  we  have 

Corollary  3.1. 

Tn(l^,cr)>Tn(iJ.,a),  n=  1,2,  ...,nsc.  (37) 

□ 

We  define  the  residual  ?■»,*;(■;  /a,  a)  :  X^0  — >  R.  for  1  <  k  <  Ki,  and  1  <  *  <  /  as 
ri,k{v;  /c,  cr)  =  -Bi(ipiik  +  bitk(n,a),v,  /x,cr),  Vi;  G  X^, 


and  the  error  bound  [13] 


|| 6j,fe(/x,cr)  -  6i,fc(M,o-)||x,i  <  A i,k(n,cr)  = 


where 


^  ,  x  r%,k{v\  H,  cr) 

sup  -W-. : - 

T,exNn  Fllx.i 


is  the  dual  norm  of  the  residual,  and  a)  is  a  lower  bound  for  the  coercivity  constant 


i(/i,cr)=  inf 


Bi(w,w;fj,,cr) 


exA,  M 


(38) 


X,i 
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Such  a  lower  bound  can  be  derived  by  hand  for  simple  parametric  dependences,  otherwise  it  can  be 
computed  using  a  successive  constraint  linear  optimization  method  [14].  It  is  also  possible  to  compute 
during  the  Offline  stage  different  coercivity  constants  of  (/n*,,  cr^)1 11  for  given  sampling  parameters  (/ 'i a*,), 
and  then  a  lower  bound  over  the  full  parameter  range  is  obtained  by  taking  the  minimum  over  the 
sampling  parameters  afB  =  min*,  aA(/Xfc,  cxfc)-  This  is  the  latter  option  that  is  used  in  the  results  section 
of  this  paper. 

We  now  assume  that  Hypothesis  2.1  holds.  Suppose  we  have  found  an,  the  nth  “shift”  such  that 
B(/x,  an)  has  a  zero  eigenvalue,  i.e.  we  have  r„(/i,  an)  =  0.  Then  our  RB-based  approximation  to  the  nth 
eigenvalue  is  A „(/x)  =  an.  We  will  now  develop  a  first  order  error  estimator  for  Tn(n,an).  We  have 

B(/x,  On)  =  r„(/x,  cr„)A(/x,  crn)Y{n,  a„), 

and  hence  with  B {n,an)  =  B(/x,  an)  +  <5B(/z,  <rn),  A(/x,  <rn)  =  A(/x,  <j„)  +  <5A(/x,  un),  V (/x,cr„)  =  V(/x,  cr„)  + 
<5V(/x,  CTn),  we  obtain 

(©(Ml  o-n)  +  <®(/x,  CT„))( V(/X,  CT„)  +  C>V(/X,  cr„))  =  Tn(/X,  crn)(A(/x,  CT„)  +  <5A(/x,  crn))(Y(/x,  CT„)  +  <5V(/x,  CT„)). 

(39) 

Expansion  of  the  above  expression  yields 

B(/x,  crn)5V(n,  an)  +  <5B(/z,  crn)V(/x,  cxn)  +  £B(/x,  an)SV(fj,,  an)  = 

Tn(ii,an)(A(fx,an)Y{fi,an)  +  A (/x,  cr„)JV(/x,  crn)  +  <5A(/z,  cr„)Y(/x,  cr„)  +  6A(/x,  crn)<5V(/x,  tr„)),  (40) 

where  the  identity  B(/x,  crn)V(/x,  crn)  =  0  has  been  employed.  We  then  multiply  through  by  V(/x,  <x„)T  and 
note  that  V(/x,  cr„)TB(^,  an)5V(n,  an)  =  (5V(^,  ct„)tB(/x,  cr„)V(/x,  crn)  =  0,  V(/x,  cr„)TA(/x,  ct„)V(/x,  an)  = 
1  and  neglect  higher  order  terms  to  obtain 

t„(/x,  Cn)  ~  V(/x,  ct„)T(5B(^,  an)Y(n,  an).  (41) 


We  then  have  the  following  bound 

I  Ki  I  Kj 

|Y(M,cr„)T<5B(^,cr„)V(^,crn)|  <  EEEE  |Yp4  (fe)  (/x,  rrn)  |  Aj  k  (/x,  cr„)  A  ji  (/x,  (J„)  p)  (/x,  <t„)  |  —  A  (/./.  cr„). 

i=l  /c=l  j=l  £=1 

(42) 

From  Proposition  2.1  part  (iii),  we  can  only  infer  eigenvalues  of  (1),(2)  when  r„(/x,  cr)  =  0,  hence  (42) 
does  not  give  us  a  direct  bound  on  the  error  of  An(/x).  However,  with  the  assumption  that  A  (/x,  <rn)  — ►  0 
in  the  limit  as  TV  — ►  oo,  we  see  that  r„(/x,  <rn)  — »  0  and  hence  asymptotically  we  have  that  An(/x)  converges 
to  A n(/x).  Moreover,  we  can  develop  an  asymptotic  error  estimator.  From  Proposition  2.2,  we  have 

rn(i tx,  An(/x))  «  r„(/x,  A„(/x))  +  (An(/x)  -  Aw(/x))  dTn^^n^ 

_  Ara(lt)  ~  A  re  (a*) 

A„(/x) 

Combining  (42)  and  (43)  gives  the  following  asymptotic  (relative)  error  estimator 


I  A»  (a*)  -^n(Af)l 
An(/x) 


< 


A(/x,<7n). 


1by  computing  the  minimum  eigenvalue  corresponding  to  (38). 
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(44) 


4.  Port  reduction 


4-1-  Empirical  mode  construction 

In  practice,  for  the  basis  functions  of  the  port  space  Z^f .  we  use  a  simple  Laplacian  eigenmode 
decomposition,  corresponding  to  the  eigenfunctions  C,Vtk  of  the  following  eigenproblem 

f  VC p,k  •  Vu  =  A  Ptk  [  Cp.fcW,  Vw  £Z*,  1  <k<  A (45) 

Jrp  Jrp 

We  can  truncate  the  Laplacian  eigenmode  expansion  in  order  to  reduce  MJ,  -  often  without  any  significant 
loss  in  accuracy  of  the  method.  However,  we  can  obtain  better  results  by  tailoring  the  port  basis  functions 
to  a  specific  class  of  problems.  A  strategy  for  the  construction  of  such  empirical  port  modes  is  presented 
in  [10].  We  briefly  describe  this  strategy  here  and  refer  the  reader  to  [10]  for  further  detail. 

A  key  observation  is  that,  in  a  system  of  components,  the  solution  on  any  given  interior  global  port 
is  “only”  influenced  by  the  parameter  dependence  of  the  two  components  that  share  this  port  and  the 
solution  on  the  non-sliared  ports  of  these  two  components.  We  shall  exploit  this  observation  to  explore 
the  solution  manifold  associated  with  a  given  port  through  a  pairwise  training  algorithm. 

Algorithm  1  Pairwise  training  (two  components  connected  at  global  port  Tp) 

<Spair  =  0- 

for  71  —  1 ,  .  .  .  ,  A samples  do 

Assign  random  parameters  a  £  [0,  <rmax]  and  /i,  £  1 ~>i  to  component  i  =  1,2 
(note  the  value  of  a  is  the  same  for  both  components). 

On  all  non-shared  ports,  assign  random  boundary  conditions. 

Solve  B(u(p,  cr),w,  p;  a)  =  0,  \/v  £  X$(p,  a) 

Extract  solution  itjpp  on  shared  port. 

Subtract  the  average  and  add  to  snapshot  set: 

Span-  S  U  ^w|rp  ~  |p  |  J  «|rp 

end  for 


To  construct  the  empirical  modes  we  first  identify  groups  of  local  ports  on  the  components  which 
may  interconnect;  the  port  spaces  for  all  ports  in  each  such  group  must  be  identical.  For  each  pair  of 
local  ports  within  each  group  (connected  to  form  a  global  port  Tp),  we  execute  Algorithm  (1):  we  sample 
this  1  =  2  component  system  many  (Asampies)  times  for  random  (typically  uniformly  or  log-uniformly 
distributed)  parameters  over  the  parameter  domain  and  for  random  boundary  conditions  on  non-shared 
ports.  For  each  sample  we  extract  the  solution  on  the  shared  port  rp;  we  then  subtract  its  average  and 
add  the  resulting  zero-mean  function  to  a  snapshot  set  S'pair.  Note  that  by  construction  all  functions  in 
Spair  are  thus  orthogonal  to  the  constant  function. 

Upon  completion  of  Algorithm  1  for  all  possible  component  connectivity  within  a  library,  we  form  a 
larger  snapshot  set  Sgroup  which  is  the  union  of  all  the  snapshot  sets  S'pair  generated  for  each  pair.  We 

12 


then  perform  a  data  compression  step:  we  invoke  proper  orthogonal  decomposition  (POD)  [15]  (with 
respect  to  the  L2(Tp)  inner  product).  The  output  from  the  POD  procedure  is  a  set  of  mutually  L2(Vp)- 
orthonormal  empirical  modes  that  have  the  additional  property  that  they  are  orthogonal  to  the  constant 
mode. 

4-2.  Port-reduced  system 

In  practice  we  use  SCRBE  -  RB  approximations  for  the  bubble  functions  -  but  as  we  will  see  in 
the  result  section,  the  error  introduced  by  RB  approximation  is  very  small  and  negligible  compared  to 
the  error  due  to  port  reduction.  As  a  consequence,  we  describe  the  port  reduction  procedure  starting 
from  the  “truth”  static  condensation  system  (30),  but  we  will  in  practice  apply  the  port  reduction  to 
the  SCRBE  system  (34).  We  recall  that  on  port  p  the  full  port  space  is  given  as 

<  =  {CP, CP,Atr}.  (46) 

For  each  port,  we  shall  choose  a  desired  port  space  dimension  ua,p  such  that  1  <  ua,p  <  A fp.  We  shall 
then  consider  the  basis  functions  Cfc>  1  <  k  <  ua,p,  as  the  active  port  modes  (hence  subscript  a);  we 
consider  the  ni)P  =  Afp  —  ua,p  remaining  basis  functions  nA,p  +  1  <  k  <  A fp,  as  inactive  (hence 
subscript  i).  Note  that  span{£Pil> . . .  ,  CP,nA,P}  C  Zp-  We  then  introduce 

nr  nr 

nA  =  Y  nA,P.  ni  =  Y  nL>  (47) 

p— 1  p— 1 

as  the  number  of  total  active  and  inactive  port  modes,  respectively;  and  n$c  =  nA  T  n\  is  the  total 
number  of  port  modes  in  the  non-reduced  system. 

Next,  we  assume  a  particular  ordering  of  the  degrees  of  freedom  in  (30):  we  first  order  the  degrees 
of  freedom  corresponding  to  the  nA  active  system  port  inodes  and  then  by  the  degrees  of  freedom 
corresponding  to  the  n\  inactive  system  port  modes.  We  may  then  interpret  (30)  as 


BAA(foo) 

BAi(a«,  o) 

V(/i,  a)  =  t{h,(j) 

AAa(m,  o) 

AAi  (n,cr) 

V(/x,  cr), 

(48) 

Bia(m,  cr) 

®h(Mj  O’) 

Aia(m,  o) 

An  (/r,  cr) 

where  the  four  blocks  in  the  matrices  correspond  to  the  various  couplings  between  active  and  inactive 
modes;  note  that  Baa(m)  G  MnAXriA  and  that  Bn(/i)  €  ]R"lXrai.  Our  port-reduced  approximation  r(/u,cr) 
shall  be  given  as  the  solution  to  the  nA  x  nA  system 

Baa(m>  ct)Va  (/X,  cr)  =  r(/u,  (j)Aaa(m,  ct)Va(m;  o'), 

Ya(m.  o)t  AAa(m,  o)Ya(m,  o)  =  1  (49) 

in  which  we  may  discard  the  (presumably  large)  Bn(/i,o)  and  An (/x,  cr)  blocks;  however  the  Bia(^;o)- 
block  is  required  later  for  residual  evaluation  in  the  context  of  a  posteriori  error  estimation. 
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4-3.  Port  reduction  error  estimator 

We  put  a  ~  on  top  of  all  the  port  reduced  quantities.  Suppose  we  have  found  an  such  that  Tn{fi ,  an)  = 
0  with  eigenvector  of  size  nsc  hr  the  non-reduced  space 

VA  ,n{P,Vn) 

0 

We  can  expand  V„(/q  an)  in  terms  of  the  eigenvectors  Vm(/x,  an)  of  the  non  reduced  space 

nsc 

v„(/r,cr„)  =  ^2  am(p,an)Vm(p,an). 

m— 1 

Since  Tn(/j,,an)  =  0,  we  can  reasonably  assume  that  |rn(/i,  an)\  =  min  crn)\.  We  now  look  at 

l<m<nsc 

the  following  residual 

nsc 

lB(/i,  CTn)Vn(/i,  (Jn)  —  ^  ^  ^"n)^m(/^5  ®n) 

m=  1 
nsc 

—  ^  ^  (M?  &n)l~m  (Ah  (M?  ^n) 

m=l 

so  using  the  A (p,an)  orthogonality  of  the  Vm(/q<rn)  we  obtain 

^SC 

||®(M!  (Jra)Yn(/-i,  Cr„)  |Ia(/j,ct?1)-i  =  ^  '  |'rrra(At)  ^ra)  I  1 1  C^m  (A4?  <7ra)Yrra(^i,  °Vi)  || 

m— 1 

nsc 

—  l'rra(Al)  °Vi)|  ^  '  ||am(/r,  cr„)A(/z,  an)Ym{/a,  (T„)\\Hllt<7n)-i 

m—1 

nsc 

=  |r”n(Al5C7ra)|  II  ^  ]  am{Pi  Cra)A(/Z,  <7n)Vm(/i,  0'n)||A(pi(Tji)-i 
m—1 

=  |T„(Ai,o-„)|2. 


V„(/i,  tT„)  = 


where  we  use  the  Euclidean  norm  derived  from  the  A(^,  cr„)  1  scalar  product.  We  thus  obtain  the 
following  error  bound 


A(/i, 


^"n)  ||  A(/x,<t71)-1  ^  ^n)  |  • 


Finally,  we  recover  an  error  estimator  for  the  eigenvalue  Xn(fJ.)  of  the  original  eigenproblem.  Assuming 
A„(/j)  is  close  to  A n(/z),  we  can  then  use  Proposition  2.2  as  in  (43),  and  we  get  the  relative  error  estimator 


I A  n{tl)  ^(aOI 

A  n(fJ.) 


<  A {n,an). 


It  is  important  to  note  that  A(/i,  an)  will  only  decrease  linearly  in  the  residual,  whereas  the  actual 
eigenvalue  error  is  expected  to  decrease  quadratically  in  the  residual.  This  is  due  to  the  fact  that  port 
reduction  can  be  viewed  as  a  Galerkin  approximation  over  a  subspace  of  the  skeleton  space  Xg(/a,  a),  and 
in  that  framework  several  a  priori  and  a  posteriori  error  results  demonstrate  the  quadratic  convergence 
of  the  eigenvalue  [16].  As  a  consequence  the  effectivity  of  the  error  estimator  A(p,,an)  is  expected  to 
decrease  with  riA.p- 
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Note  that 


M(p,an)Yn(p,crn)  = 


0 


®Ia(/^?  ^"n)^A,n(/^5  &n)  J 

and  so  the  computation  of  the  residual  requires  the  additional  assembly  of  Bia(m,  crn),  which  does  not 
generate  an  important  extra  computation  since  in  practice  we  will  consider  ua  n\.  On  the  contrary, 
the  computation  of  the  norm  ||  •  Ha^oa)-1  requires  the  assembly  and  inversion  of  A(p,an),  the  full 
Schur  complement  stiffness  matrix,  which  would  potentially  eliminate  any  speed-up  obtained  by  the  port 
reduction.  This  computational  issue  is  resolved  by  using  an  upper  bound  for  ||  •  ||  ao.oa)-1  which  is  based  on 
a  non-conforming  version  A'  (p.  an)  of  the  stiffness  operator  and  a  parameter  independent  preconditioner: 
the  former  permits  online  computation  of  small  matrice  inverses  locally  on  each  component,  and  the  latter 
allows  us  to  precompute  non-reduced  matrices  and  their  Cholesky  decompositions  in  an  offline  stage. 
The  entire  procedure  is  described  in  detail  in  [10]. 


5.  Numerical  results 


5.1.  Linear  elasticity 

We  consider  linear  elasticity  in  a  non-dimensional  form:  we  nondimensionalize  space  with  respect  to 
a  length  do  which  will  correspond  to  the  beam  width  in  the  following,  we  nondimensionalize  the  Young’s 
modulus  E  with  respect  to  a  reference  value  Eq,  and  we  nondimensionalize  time  with  respect  to  y J 
where  p  is  the  mass  density.  The  non  dimensional  linear  elasticity  free  vibration  equation  then  reads 

d2U 


- AU  = 


dt 2  ’ 


(50) 


where  A  is  a  linear  second  order  differential  operator  in  space  and  U ( x ,  t)  is  the  displacement  vector. 
Assuming  that  the  free  vibration  solution  is  of  the  form  U(x,t )  =  u(x)cos(uit) ,  the  problem  is  equivalent 
to  solving  the  eigenproblem 

Au  =  uj2u.  (51) 


In  variational  form,  the  operator  A  corresponds  to  the  bilinear  form  [13] 

a(w,v;p)  =  /  Cijki (p)eij (w)eki (v)  (52) 

Jn 

where  we  assume  summation  on  repeated  indices;  a(-,  •)  is  defined  on  the  space  of  admissible  displacements 

V  =  {v  =  (vi,V2,Vo)\vi  €  H1(Ll)-,Vi  =  0  on  T0  C  dLl},  and  eij(v)  =  ^(i diVj  +  djVi).  We  will  consider 

isotropic  materials,  in  which  case  the  coefficents  Cijki(p)  are  functions  of  only  two  parameters:  Poisson’s 

ratio  v  and  Young’s  modulus  E.  In  the  following  we  always  fix  v  =■  0.3,  and  allow  E  to  vary,  hence  E  is 

part  of  the  vector  of  parameters  p.  More  precisely,  the  parametric  dependence  reads 

Eta  E 

CijklM  =  (1  _|_  __  2V)  +  2(1  _|_  p\ 

We  also  define  the  mass  bilinear  form 


m(w,  v ;  p)  = 

15 


(53) 


Figure  1:  The  component  library:  a  beam  (left)  and  a  connector  (right).  Components  can  connect  at 
square  ports  shown  in  red.  All  ports  have  the  same  shape  and  same  discretization. 


The  eigenproblem  in  variational  form  finally  reads:  find  A(/r)  £  R>o  and  u(/r)  £  V  such  that 

a(u(p),v,  fj)  =  A (/j,)m(u(n),v\  n)>  Vu  £  V,  (54) 

m(u(n),u(ii);n)  =  1-  (55) 

Note  that  A(/z)  =  w2(/r)  -  the  eigenvalue  is  the  frequency  squared. 

5.2.  Component  library 

We  consider  a  linear  elasticity  library  of  two  components  shown  in  Figure  1:  a  beam  and  a  connector. 
The  FE  hexahedral  meshes  are  shown  in  Figure  1,  and  in  all  the  following  we  use  first  order  approximation 
with  trilinear  elements.  The  components  can  connect  at  square  ports  of  dimension  lxl  with  A = 
3  x  36  =  108  degrees  of  freedom.  The  beam  has  two  parameters:  the  Young’s  modulus  E  £  [0.5,  2]  and 
the  length  scaling  s  £  [0.5,2],  where  the  beam  is  of  length  5s.  The  connector  has  one  parameter,  the 
Young’s  modulus  E  £  [0.5,  2].  Finally,  for  the  shift  parameter  a,  we  consider  the  range  [0,  0.01],  based  on 
the  fact  that  the  local  minimum  eigenvalues  of  the  two  components  are  larger  than  0.01  for  the  previous 
E  and  s  parameter  ranges.  For  each  component,  we  build  RB  bubble  spaces  of  size  N  =  10  using  a 
Greedy  algorithm  [17],  for  the  parameter  ranges  previously  defined.  We  also  perform  a  pairwise  training 
for  the  component  pair  beam-connector  to  build  empirical  port  modes  as  described  in  Section  4.1;  and 
we  build  a  parameter  independent  preconditioner  (necessary  for  the  computation  of  A)  using  parameter 
values  E  =  0.5  and  s  =  0.5. 

5.3.  Simple  beam 

We  first  present  a  simple  example  where  we  compare  with  beam  theory  to  demonstrate  that  the 
FE  resolution  is  adequate  and  that  we  capture  the  different  modes.  We  consider  a  clamped-clamped 
uniform  beam  of  square  section,  with  thickness  d  =  1  and  length  L  =  40,  and  Young’s  modulus  E  =  1. 
Table  1  presents  the  first  eight  eigenvalues  obtained  by  different  methods:  Euler  Bernoulli  model  [18], 
Timoshenko  model  [18],  global  FEM  and  SCRBE  with  and  without  port  reduction  (in  which  the  beam 
is  constructed  as  the  concatenation  of  eight  beam  components).  The  eigenvalues  (which  we  recall  are 
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Figure  2:  The  fourth  (pure  bending)  and  fifth  (pure  torsion)  beam  eigenmodes. 


the  frequencies  squared)  are  quite  small  as  the  beam  is  of  large  aspect  ratio.  The  SCRBE  results  are 
obtained  by  connecting  eight  beam  components  together  with  length  parameter  s  =  1,  using  RB  spaces 
of  size  N  =  10;  no  port  reduction  corresponds  to  ua,p  =  A =  108,  and  for  port  reduction  we  use 
riA.p  =  20  active  port  modes.  The  global  FEM  results  are  obtained  using  a  global  mesh  corresponding  to 
eight  beam  component  meshes  stitched  together,  hence  SCRBE  and  FEM  are  based  on  the  same  mesh 
and  FE  resolution. 

We  observe  that  the  beam  models  do  not  capture  some  eigenvalues;  these  correspond  to  torsional 
modes  that  are  not  taken  into  account  in  Euler  Bernoulli  and  Timoshenko  models  which  consider  only 
bending  displacement.  Note  that  for  a  beam  with  a  square  section,  the  bending  and  torsion  is  decoupled 
and  the  eigenmodes  are  either  pure  bending  or  pure  torsion  (see  Figure  2).  For  the  modes  that  are  pure 
bending  (Ai;  A2,  A3,  A4,  A6,  A7),  we  observe  a  good  agreement  between  all  methods.  Note  that  it  is  well 
known  that  Euler  Bernoulli  is  better  for  long  wavelength  and/or  slender  beams;  Timoshenko  is  better 
for  shorter  wavelength  and/or  shorter  beams.  Not  surprisingly,  the  FE  (and  SCRBE)  eigenvalues  are 
closer  to  Euler  Bernoulli  for  lower  modes  and  closer  to  Timoshenko  for  higher  modes.  The  SCRBE  (with 
or  without  port  reduction)  and  global  FEM  give  results  that  have  an  actual  relative  difference  less  than 
10~4.  For  the  SCRBE  without  port  reduction,  we  also  give  the  relative  error  estimate  A  in  Table  1, 
which  corresponds  to  the  relative  error  between  the  SCRBE  and  the  “truth”  static  condensation:  it  is 
at  most  10~6,  which  confirms  that  the  error  introduced  by  RB  is  negligible.  For  the  SCRBE  with  port 
reduction  and  ua.p  =  20,  the  relative  error  estimate  is  A  and  corresponds  to  the  relative  error  between 
SCRBE  with  and  without  port  reduction:  it  is  about  10-2,  which  overstimates  the  actual  relative  error, 
but  nonetheless  indicates  a  very  good  agreement  between  SCRBE  eigenvalues  with  and  without  port 
reduction.  We  also  observe  that  the  SCRBE  does  capture  all  the  torsional  modes.  Note  finally  that 
the  SCRBE  eigenvalues  are  obtained  using  a  root  finding  algorithm:  in  practice  we  set  the  tolerance2 
to  10~10  as  this  is  a  couple  orders  of  magnitude  smaller  than  the  RB  relative  error  estimator  A,  thus 
making  the  root  finding  error  negligible  with  respect  to  RB  error  (and  also  port  reduction  error). 

5-4-  Bridge  structure 

We  are  now  ready  to  consider  larger  systems  with  more  complicated  connections  which  will  better 
exercise  the  RB  and  port  reduction  capabilities.  Towards  this  end,  we  consider  a  system  of  30  components, 
corresponding  to  a  bridge  structure,  shown  in  Figure  3  with  its  first  three  eigenmodes;  in  that  case  we 
choose  E  =  0.5  and  s  =  1  for  all  components.  In  the  following,  we  will  provide  systematic  analysis  of 

2the  tolerance  applies  to  rn (//.,  er)  =  0,  hence  it  corresponds  to  a  relative  tolerance  for  Arf(//). 
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Ai 

A2 

^3 

A4 

A5 

A6 

A7 

As 

Euler  Bernoulli 

1.6294e-05 

1.2381e-04 

4.7583e-04 

1.3003e-03 

— 

2.9016e-03 

5.6603e-03 

— 

Timoshenko 

1.6204e-05 

1.2224e-04 

4.6524e-04 

1.2560e-03 

— 

2.7622e-03 

5.2991e-03 

— 

Global  FEM 

1.6612e-05 

1.2489e-04 

4.7327e-04 

1.2708e-03 

2.0732e-03 

2.7775e-03 

5.2916e-03 

6.1912e-03 

SCRBE  nAiP  =  108 

1.6612e-05 

1.2489e-04 

4.7327e-04 

1.2708e-03 

2.0732e-03 

2.7775e-03 

5.2916e-03 

6.1912e-03 

SCRBE  nA,P  =  20 

1.6612e-05 

1.2489e-04 

4.7327e-04 

1.2708e-03 

2.0732e-03 

2.7775e-03 

5.2916e-03 

6.1912e-03 

A 

1.4418e-06 

2.0695e-07 

7.9612e-08 

9.6913e-08 

5.4576e-09 

4.1418e-07 

1.0262e-06 

8.8249e-09 

A 

5.5488e-03 

7.3845e-03 

8.4207e-03 

7.4811e-03 

3.3180e-02 

8.3262e-03 

8.9995e-03 

4.7761e-03 

Table  1:  Eigenvalues  for  a  clamped-clamped  uniform  beam  of  square  section,  with  thickness  d  =  1  and 
length  L  =  40.  The  estimators  A  and  A  correspond  relative  errors. 


Figure  3:  A  bridge  structure:  all  of  the  “open  ports”  of  beam  components  are  clamped.  The  first  three 
displacement  eigenmodes  are  shown. 
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Figure  4:  First  eigenvalue  convergence  with  respect  to  the  size  N  of  RB  spaces. 

the  RB  and  port  reduction  convergence  and  also  performance  of  the  a  posteriori  error  estimates;  finally, 
we  will  report  detailed  timings  to  confirm  the  computational  advantage  for  large  problems. 

We  first  show  in  figure  4  the  convergence  of  the  first  eigenvalue  with  respect  to  the  size  N  of  the  RB 
spaces  used  for  bubble  approximations.  Note  that  we  did  not  compute  the  eigenvalue  with  the  “truth” 
static  condensation,  because  it  would  be  very  computationally  intensive,  hence  the  reference  value  for 
A  is  the  value  obtained  with  a  global  FEM,  denoted  A fe3-  We  observe  that  we  obtain  exponential 
convergence,  hence  we  provide  a  significant  improvement  compared  to  standard  CMS  approaches.  We 
also  observe  that  the  RB  relative  error  estimator  A  is  accurate  -  it  overestimates  the  actual  error  by  at 
most  one  order  of  magnitude;  moreover,  for  a  sufficiently  large  TV,  the  RB  relative  error  estimator  A  is 
very  small,  hence  justifying  the  fact  that  we  can  neglect  the  error  due  to  RB  error  approximation  when 
introducing  port  reduction. 

We  now  fix  N  =  10  and  consider  port  reduction.  We  show  in  Figure  5  the  convergence  of  the  first 
eigenvalue  with  respect  to  the  number  of  port  modes,  for  both  the  regular  Laplacian  modes  and  the 
empirical  modes:  the  advantage  obtained  with  the  empirical  modes  is  obvious,  and  we  also  observe  that 
A  does  not  converge  as  fast  as  the  true  error,  which  is  due  to  the  fact  that  it  is  only  a  linear  function  of 
the  norm  of  the  residual,  as  explained  in  Section  4.3. 

We  now  consider  computational  savings  compared  to  a  global  FEM  approach.  We  report  in  table  2  the 
first  eigenvalue  computation  time  and  error  bounds  with  and  without  port  reduction;  in  table  3  the  same 
results  are  reported  for  the  tenth  eigenvalue  computation.  Note  that  we  use  the  Krylov-Schur  method 


’'Tlie  global  FEM  eigenvalue  is  not  expected  to  be  exactly  the  same  as  what  would  be  obtained  with  FE  static  con¬ 
densation  -  in  theory  they  should  be  the  same,  but  the  different  computational  paths  lead  to  different  numerical  results  - 
hence  it  explains  why  ^  does  not  converge  to  zero,  and  why  A  gets  smaller  than  ^ ^  ^  for  N  big  enough. 
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Figure  5:  First  eigenvalue  convergence  with  respect  to  the  number  71a.p  of  active  port  modes: 
responds  to  Laplacian  eigenmodes  (“Lap”),  red  corresponds  to  empirical  modes  (“Emp”).  The 
eigenvalue  A  is  the  one  obtained  for  riA.p  =  =  108. 

with  the  shift  and  invert  transform,  provided  by  the  SLEPC  library  [19].  For  SCRBE,  we  use  a  direct  LU 
solver  for  matrix  inversion  (thanks  to  the  moderate  number  of  dofs);  for  FE  we  use  an  iterative  solver: 
conjugate  gradient  with  incomplete  Cholesky  preconditioner  (1  cpu)  or  Jacobi  preconditioner  (16  cpu). 
For  the  first  eigenvalue,  without  port  reduction,  the  SCRBE  yields  a  speed-up  factor  of  60  compared 
to  FE  with  1  cpu,  and  a  speed-up  factor  of  11  compared  to  FE  with  16  cpu.  For  the  tenth  eigenvalue, 
the  decrease  in  computational  time  is  even  larger  -  112  compared  to  FE  with  1  cpu  and  23  compared 
to  FE  with  16  cpu.  Computing  eigenvalues  that  are  not  at  an  extremity  of  the  spectrum  is  a  more 
difficult  problem,  and  in  that  case  the  computational  advantage  of  the  SCRBE  is  even  more  obvious.  We 
observe  that  the  relative  difference  between  the  SCRBE  and  FE  eigenvalues  is  less  than  10-6.  We  now 
add  port  reduction:  we  retain  only  some  of  the  port  (empirical)  modes  for  the  eigenproblem  solution. 
When  retaining  about  a  fifth  of  the  port  modes  («a,p  =  20),  the  speed-up  factors  for  the  first  and  tenth 
eigenvalues  are  960  and  2500  compared  to  FE  with  1  cpu,  180  and  520  compared  to  FE  with  16  cpu. 
The  SCRBE  eigenvalues  are  again  within  a  10”6  relative  distance  from  the  FE  eigenvalues,  and  the  port 
reduction  relative  error  estimator  predicts  a  1%  relative  error  for  the  port-reduced  eigenvalues  ■  despite 
overestimation  of  the  true  error,  this  1%  relative  error  estimate  is  more  than  sufficient  in  an  engineering 
context.  Note  that  the  SCRBE  can  be  easily  parallelized  thanks  to  the  component  assembly  description, 
but  in  this  paper  the  SCRBE  computations  are  done  on  a  single  CPU,  hence  the  comparison  of  global 
FEM  with  16  cpu  to  SCRBE  with  1  cpu  is  strongly  biased  in  favor  of  FEM. 

For  our  formulation  to  remain  well  posed,  we  can  only  capture  eigenvalues  that  lie  in  the  a  range 
[0,0.01].  Despite  this  limitation,  we  are  able  to  compute  up  to  A50  =  0.00994974,  which  is  more  than 


blue  cor- 
reference 
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number  dofs 

1st  eigenvalue 

A  {/I,  a) 

A  (h,ct) 

timing 

Global  FEM  (1  cpu) 

AT  =  267792 

0.000108903 

32  min 

Global  FEM  (16  cpu) 

A f  =  267792 

0.000108903 

6  min 

SCRBE  tia.p  =  Afp  =  108 

nA  =  nsc  =  6480 

0.000108903 

1.36e-08 

32s 

SCRBE  tia.p  =  20 

nA  =  1200 

0.000108903 

1.16% 

2s 

Table  2:  First  eigenvalue  computation  for  a  bridge  assembled  as  30  components  (Figure  3).  E  =  0.5  and 
s  =  1  everywhere. 


number  dofs 

10th  eigenvalue 

A (At,  cr) 

A(/i)  c) 

timing 

Global  FEM  (1  cpu) 

AT  =  267792 

0.00194414 

125  min 

Global  FEM  (16  cpu) 

A f  =  267792 

0.00194414 

26  min 

SCRBE  tia.p  =  Afp  =  108 

n  a  =  nsc  =  6480 

0.00194414 

2.50e-08 

67s 

SCRBE  nA,P  =  20 

nA  =  1200 

0.00194415 

1.67% 

3s 

Table  3:  Tenth  eigenvalue  computation  for  a  bridge  assembled  as  30  components  (Figure  3).  E  =  0.5 
and  s  =  1  everywhere. 

enough  in  this  structural  analysis  case  where  only  the  lowest  frequency  modes  are  likely  to  resonate. 

We  demonstrated  that  SCRBE  has  a  computational  advantage  relative  to  FE,  but  the  same  decrease 
in  computational  time  could  in  theory  be  obtained  with  CMS.  One  crucial  advantage  of  SCRBE  with 
respect  to  CMS  (in  addition  to  convergence  rate)  is  its  flexibility  with  respect  to  parameter  variations. 
Thanks  to  the  RB  approximations  at  the  component  level,  we  can  modify  the  component  parameters  and 
recompute  the  eigenproblem  solution  seamlessly.  We  show  in  Figure  6  the  third  eigenmode  for  different 
parameter  variations:  we  can  modify  some  of  the  beam  lengths  (Figure  6a),  or  we  can  make  one  half  of 
the  bridge  stiffer  than  the  other  (Figure  6b). 


Figure  6:  Comparison  of  the  third  eigenmode  for  different  parameters.  Left:  for  the  middle  beams, 
s  =  0.7;  for  the  beams  adjacent  to  the  middle  beams,  s  =  1.3;  for  all  other  beams  and  all  support  beams, 
s  =  1.  Right:  in  the  first  half  of  the  bridge  E  =  2,  in  the  second  half  E  =  0.5. 
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Figure  7:  A  notched  beam  component. 


Original  system 

Notched  beam 

Disconnected  beam 

Ai 

0.000108903 

0.000109415 

0.000107512 

A2 

0.000192330 

0.000192394 

0.000152107 

A3 

0.000351258 

0.000351093 

0.000336273 

a4 

0.000502429 

0.000505191 

0.000365281 

A5 

0.001227793 

0.001206475 

0.000874600 

Table  4:  Eigenvalue  comparison  when  introducing  defect  in  a  beam. 

Finally,  one  could  argue  that  the  examples  presented  in  this  section  have  relatively  simple  geometries 
and  that  the  same  results  could  be  obtained  by  assembling  components  based  on  beam  models,  without 
solving  the  full  linear  elasticity  equations.  But  one  of  the  advantages  of  our  approach  (and  the  CMS  in 
general)  is  its  ability  to  handle  arbitrary  component  shapes,  allowing  to  tackle  problems  that  are  out  of 
reach  for  beam  and/or  plate  models.  As  an  example,  we  introduce  a  notched  beam  component  (Figure  7) 
in  our  original  bridge  system:  consequently,  a  new  mode  appears  between  the  ninth  and  tenth  eigenmodes 
of  the  original  system,  as  seen  in  Figure  8.  This  change  of  behaviour  would  be  impossible  to  predict  only 
using  beam  models.  To  emphasize  the  latter,  we  compare  our  notched  beam  case  to  a  “disconnected 
beam”  case  (with  zero-thickness  crack,  Figure  9c)  that  would  be  the  only  possible  modeling  of  a  defect 
with  beam  models:  table  4  clearly  shows  that  a  notched  beam  introduce  subtle  changes  in  the  spectrum 
whereas  a  disconnected  beam  completely  modifies  the  spectrum. 
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